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Abstract 

We study the emergence of the collective spatio-temporal macro- 
scopic properties of the immune system, by representing individually 
the elementary interactions between its microscopic components (anti- 
bodies, antigens, cytokines). The results of this detailed explicit anal- 
ysis are compared with the traditional procedure of averaging over 
individuals and representing the collective dynamics in terms of den- 
sities that obey partial differential equations (PDE). The simulations 
show even for very simple elementary reactions the spontaneous emer- 
gence of localized complex structures. In turn the effective dynamics 
of these structures affects the average behavior of the system in a very 
decisive way: systems which would according to the differential equa- 
tions approximation die, display in reality a very lively behavior. Our 
conclusions are supported both by explicit microscopic simulations 
and by analytic calculations. As the optimal method we propose a 
mixture of microscopic simulation (MS) systems describing each reac- 
tion separately, and average methods describing the average behavior 
of the agents. 

1 Introduction 

The immune system has become lately a classical example of emergence 
It is an archetype for biological systems whose behavior cannot be directly 
deduced from their molecular and cellular properties 0. 

As in many other biological systems, a large amount of microscopic data 
has been gathered. This data by itself is insufficient to explain the macro- 
scopic behavior of the immune system. 

Like any complex system, the immune one, requires classically a few com- 
plementary ways for describing it. 

• The micro-biological level - molecules and cells. 

• The generic properties that give the broad picture of what makes this 
system tick. For instance the "clonal selection": the random genetic 
changes of the antibodies and their selection by the affinity with ap- 
propriate antigens. || 

• The macroscopic mechanistic level. For example the aggregated ac- 
tion of groups of cell and cytokine: TH1 vs TH2, idiotypic vs anti 
idiotypic... 
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• The complex network of interactions. This represents the emergence 
of the immune response in terms of a highly connected (e.g neural) 
network. In particular, one uses neural networks to implement this 
view 

• The reaction-diffusion (partial) differential equations level that describes 
the IS ClS db nonlinear chemical system 0. 

We are proposing here a methodology, which has the capacity to express 
much of the successes of the other methods while transcending their lim- 
itations and in particular the limitations which prevented the appropriate 
understanding of the central issue: the emergence of complex deterministic 
macroscopic functions out of simple stochastic microscopic interactions . 

The methodology is based on recognizing the discreetness and spatial in- 
homogeneity of the biological units || |9| and tailoring a specific model for 
each biological situation by representing the objects and interactions appro- 
priate for each scale. 



More precisely, we introduce microscopic simulation (MS) (T^, [II], |I~2| 
methods, which on one hand represent directly the basic elements of the 
immune system but on the other hand allow the identification of the effective 
macroscopic objects relevant to the collective dynamics of the system. 

In our MS we enact in the computer each reaction separately. More 
precisely we compute at each lattice point the probability of each reaction 
and then perform reactions randomly according to these probabilities. 

We first show that very simple microscopic interactions may lead system- 
atically (when studied with MS) to highly nontrivial effects which one would 
not have guessed using any of the classical methods. 

We describe the resulting new phenomena and show how they are relevant 
for various biological applications (e.g. the generation of complex spatio- 
temporal localized self-organized adaptive structures). 

In other words we show that systems that would be judged too simple to 
constitute an adaptive network and too linear to present non-homogeneous 
solutions in the differential equations context, do in fact emerge well de- 
fined/individualized macroscopic objects, which can move around, act in 
complex ways, and display adaptive behavior. 

We present here our Microscopic Simulation and hybrid (discrete-continuum) 
models as a generic technology for representing in a direct intuitive way the 
microscopic knowledge on the immune system and for extracting non-trivial 
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macroscopic information out of them (including the emergence of adaptive 
behavior and self-organization). 

Our results are at variance with the predictions of the classical modeling 



methods ||15|| , but are in encouraging agreement with many known IS features 
which were long waiting for an explanation. 

In chap 2 we study the macroscopic effect of the proliferation of cells 
induced by a macroscopically homogeneous, but microscopically inhomoge- 
neous signal distribution. We show that the proliferation rate predicted by 
the MS is many orders of magnitude larger than the one expected from ODE. 
In particular one obtains proliferation in conditions in which the naive con- 
tinuum limit predicts decay. 

In chap 3 we identify the source of the difference between MS and ODE 
and propose hybrid models (HM) that bridge the gap between the 2 limits. 

In chap 4 we introduce additional reactions, which limit the growth of the 
proliferating entities. We show that in certain conditions the macroscopic 
inhomogeneity is even increased by the effect of competition. 

Most of the effects analyzed in the four first chapters are due to the 
combination of auto-catalysis and of fluctuations. 

In chap 5 We study the effects of discreetization in the presence of negative 
feedback loops. 

We find that some effects predicted by the continuum simulations |TB|, [T7| 
(but not so frequent in nature) are not present when the system is analyzed 
using MS (e.g. the spurious oscillatory behavior of certain regulatory pairs). 

In conclusion we are opening the way to a wide front re-evaluation of the 
dynamics that is supposed to emerge out of the known microscopic elemen- 
tary interactions of the immune system. We provide the main effects to be 
re-scrutinized and the particular methods to do so. Paradoxically, the new 
methods that we propose are closer to the intuition of the practitioners in the 
immunological field. Indeed the objects which MS manipulates are the very 
same they encounter in their experimental work: cells, cytokines, epitopes 
(rather than abstract approximations such as differential equations, neural 
networks etc). 



2 Dynamics of Discrete Proliferative Agents 

Proliferating entities can be described as agents A that duplicate whenever 
they meets a stimulating signal S. These agents have in general also a typical 
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death rate d. Examples of such agents are : in immunology - immune cells 
reacting to an antigen |18 ; In ecology - animal proliferating whenever they 
find food... 

We will deal with the dynamics of the stimulating signal in the next 
sections. In this section we will assume the stimulating agent S is fixed in 
time. 

This system can be described as a 2 reactions system (Figure 1). The 
first reaction is that whenever am A agent meets a stimulating signal S it 
multiplies with a probability d±. We will note it as : 

• A + S -> A + A + S. di 

The second mechanism is the death of an agent with a probability d 2 

• A -> $. d 2 

In the homogeneous limit the growth rate of A is proportional to number 
of A and S pairs. Thus the number of new A agents produced in a time 
interval At is: 

AA = SA * d 1 AT (1) 
The number of A deaths is proportional to the total A population: 

AA = -d 2 * AAt (2) 

Thus the total change in A is: 

AA = (S * di - d 2 ) * AAt (3) 

This equation becomes an ODE [19j when we take the limit At —>■ 

A = (S * d 1 - dj,) * A (4) 

The solution of the ODE x = ax is x = e at 111 911. In this case the solution is : 



A = e^ 1 -^* 1 (5) 
This solution has 3 regions of behavior (Figure 2): 

a) If the average growth rate is higher than the death rate i.e. (S * d\ — 
d 2 ) > then the A concentration grows with an exponential growth rate 
of (S * d\ — d 2 ) This exponential growth will end when other mechanisms, 
discussed in the next section, become important. 
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b) If the average growth rate is lower than the decay rate i.e. (< S > 
*di — d 2 ) < the A concentration decays to 0. 

c) If the average growth rate is equal to the decay rate i.e. (< S > 
*di — d 2 ) = the A concentration is constant in time. 

This solution does not take into account the non-homogeneous distribu- 
tion of S. In a real system the proliferation rate changes from one point 
to the other. Regions with a high S concentration have a high prolifera- 
tion rate, while regions with a low S concentration have a low (or negative) 
proliferation rate. 

In the ODE we assumed that the effective proliferation rate is equal to 
the average proliferation rate (< S > *d\ — d 2 ). This assumption can be 
very mistaken in many cases. Indeed if the interaction between S and A are 
local then the local dynamics is described by the local population A, at each 
location i. 

K = (Si * d 1 - d 2 ) * Ai (6) 
The correct averaging of equation |6| is: 

A = d x < SiAi > -d 2 A (7) 

Eq. ^transforms into Eq. f| if we can separate of < SiAi > into < Si >< 
Ai >. Thus in order for Eq. 4 to hold S and A have to be independent 
variables p0| . Eq. 4 is the very basis of any ODE approach. Yet it is 
obviously wrong in the case we present here since Ai is high precisely in the 
sites where Si is high. Thus S and A are in fact strongly correlated. 

Following this argument we expect to see differences between the micro- 
scopic representation of this system and its description by ODE. 

Indeed if the signaling agent (S) is a random discrete agent its distri- 
bution is uniform over large scales but locally its density will contain some 
granularity. The difference between the large scale 5* uniformity and its mi- 
croscopic granularity is the origin of the macroscopic differences between the 
ODE and the MS results. 

For example if (< S > *d\ — d 2 ) < (the average birth rate is lower 
than the average death rate) the ODE would yield a decaying solution but 
the MS yields an exponential growth (Figure 3). To understand this growth 
note that each point in space contains a different concentration of signaling 
agent (S). For most of the points i one will have (Si * di — d 2 ) < and 
the local population will decay to zero. However for some points one may 
have (Si * di — d 2 ) > and the local population at these points will rise with 
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an exponential rate. After a short time the contribution of the points with 
decaying population will be negligible and only the lattice points with an 
exponential growth contribute to the average population. In particular the 
population growth is dominated by the site with the highest S number. 

The previous description does not take into account the diffusion of the 
A's. The A diffusion will lead to the creation of large scale islands of high 
A population around each maxima of The 5* density. These islands will keep 
growing until they fill all space. The diffusion of the A does not change 
qualitatively the average population growth mentioned above for a very wide 
range of parameters [12| . In order to avoid complex formulae we did not rep- 
resent explicitly the diffusion term /j,AA. In all the simulations and analytic 
computations this term has been dully taken into account. We mention it in 
the sequel only when it has a major qualitative effect on the results. 



3 How Well do Different Methods Deal with 
Discreetness 

The MS helped us to uncover 2 basic features: 

• The population increase in the MS is much higher than the one com- 
puted in the homogeneous case: The ODE is exponentially decaying 
while the MS is exponentially growing. 

• Diffusion by itself is not enough to smooth A concentration distribution. 

The effect of the non-homogeneity can be observed in every experiment 
with a petri dish in which the proliferation of cells creates first local patterns 
of very high concentration even if the original distribution is homogeneous. 
The irrelevance of diffusion can be observed in many ecological systems, in 
which although the diffusion rate of each animal is much higher than its 
reproduction rate localization can still be observed. 

3.1 Microscopic Simulation vs. PDE 

One might think that using partial (rather than ordinary) differential equa- 
tions will appropriately take into account the non-homogeneity in the a dis- 
tribution, this is actually not the case: the fact that A can diffuse will 
only increase the uniformity in the A distribution, while the microscopic 
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inhomogeneity of S is lost upon expressing the S distribution as a continu- 
ous function (recall that S was spread with uniform probability across the 
space). Consequently there is no reason to expect that the partial differential 
equations acting on smooth distributions will do better than ODE. This is 
confirmed by rigorous analytical treatment. 

In the next chapter we show that allowing S to diffuse does not solve this 
problem either: the diffusion will only interchange one microscopically inho- 
mogeneous configuration with another, equally inhomogeneous one. Their 
difference as well as their very inhomogeneity is lost once one expresses the 
distribution of discrete S in terms of a continuous density function. 

3.2 Discreetization Implies Fluctuations 

The crucial role of microscopic discreetization in the emergence of the macro- 
scopic inhomogeneity might be unfamiliar and even strange for the reader. 
After all, the microscopic corrections to the differential equations should 
morally lead to microscopic corrections in their solutions. That this is how- 
ever not the case is known in the theory of electric conductivity since the 
late fifties. Indeed, by solving the Schroedinger equations for periodic poten- 
tials one obtains electronic wave functions that spread over the entire space 
(which implies electric conductivity) |13[. However, at least in 2 dimensions, 
the slightest random uniform noise leads to spatially localized wave functions 
(insulator) ||14| . Similarly here the presence of microscopic discreetization of 
S induces intrinsically microscopic spatial inhomogeneity in the S distribu- 
tion. The partial differential equation (PDE) for A share the property that 
(at least in 2 dimensions) the slightest spatial inhomogeneity (of S) in the 
equations leads to macroscopic localization. We have checked that indeed the 
effect of S discreetness can be fully taken into account by working with un- 
quantized (real valued) S"s with same distribution function as the quantized 
5"s (Figure 4). continuous 

3.3 Hybrid Models 

The previous section showed that the difference between the MS and ODE 
results are only due to the microscopic non homogeneity of the S agents. 
We conceived a class of models containing elements from both ODE and 
MS. We call such models with the collective name of hybrid models (HM). 
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The particular combination of ODE and MS techniques has to be tailored 
specifically for each biological situation. 

In the present case we compute the local A concentration changes in a 
continuous way, but use the S distribution produced by the MS. However the 
spatial distribution of the AS is non-uniform, and in fact we don not assume 
that the A's spatial distribution is . 

In other words we replace the stochastic dynamics of the MS by one dif- 
ferential equation per lattice site, and we preserve thereby the inhomogeneity 
in the agents concentrations. 

Other system require other HM that fit their specific aspects. For example 
a system containing regions with high interaction rates, and high concentra- 
tions, and regions with low interaction rates can be modeled using a HM 
that computes the interaction rate stochastically for low concentrations, and 
deterministically for high concentrations. 

HM model has the advantage that when the number of A and S pairs 
per site is large, the local reaction rate can be computed directly rather than 
having to enact explicitly each and every A + S— > A + A + S reaction. 
Thus they have a much lower CPU cost than the MS. 

intrinsically This very simple system shows that the results one obtains 
depends crucially on the appropriate use of models, and the recognition of 
the limits of applicability of the various approximations. It draws in rough 
lines the limits between the range of applicability of each type of model. 

• ODE can be used when the diffusion rate is much higher than the 
maximum activation rate. 

• MS are precise over a very wide range but have a very high cost in 
CPU time. 

• 4D PDEs have all the severe problems related to the continuity assump- 
tion, and are not much cheaper than MS. However there is a small range 



of cases where they can provide analytical solutions [21 



• HM can reproduce most of the results produced in MS without the high 
CPU cost. They reproduce the effects of microscopic inhomogeneity. 

The HM direction can be used in an adaptive way. In particular the 
collective objects emerging form the dynamics can be treated as the effective 
elementary agents for higher (coarser scale) level of analysis. 
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4 Mechanisms Limiting the Population Growth 

The proliferation of the A agents in the models of section 2 must be limited by 
some mechanism if we do not want them to take over the universe. We have 
shown that a high death rate is not enough in order to limit the population 
proliferation. The required mechanisms can be either an external regulator 



[22], or a limit on the resources available to the proliferating cells [p3| 



4.1 Local Competition 

Many of the mechanisms limiting the A proliferation can be described as 
competition. Competition for a resource means that whenever 2 cells need 
the resource either to survive or to proliferate only one of them will get it and 
the other will die (or fail to proliferate). In ecology competition can be over 
food or water, while in the immune system competition can be over access 
to an antigen (Figure 5). 



Such a system is sometimes called a Lotka Voltera system |24| , p5| and 
can be described using a local differential equation: 

A = CL& * Ai - d 2 Ai - XAj. (8) 

The system can exist in two regimes: 

• If the proliferation rate is higher than the death rate, this system will 
proliferate until it reach the saturation level uniformly throughout the 
entire space. In such there is no difference between the MS and 
ODE results. 

• If the proliferation rate is lower than the death rate than the ODE 
would eventually converge to a uniform vanishing concentration. By 
contrast the MS will reach a steady state of islands whose population 
density is at saturation level. These islands are embedded in a back- 
ground empty of A. For high A diffusion rate the MS yields non-zero 
concentration everywhere (Figure 6) in spite of the ODE prediction. 

Paradoxically the conditions in which ODE and the MS show similar fea- 
tures are seldom realized in nature. They represent extreme cases in which a 
certain agent is capable to fill all the available space. In the animal world only 
a few species enjoy such a status (For example humans and their parasites). 
In immunology this regime corresponds to a cell type filling all the blood 
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and immune system (systemic diseases). Except for these extreme cases the 
situation in which (< S > *di — (I2) < seems to be more appropriate for 
realistic systems. 

4.2 Global Competition 

The local competition mechanism presented in the previous section assumes 
that the limiting factor that creates the competition is homogeneous on the 
same scale as the A agents concentration. (The scale of homogeneity is 
the spatial scale in which a cell type has a homogeneous concentration). 
Some reactions leading to competition have a larger homogeneity scale than 
the factor they are limiting. For example if cells proliferate and require for 
their proliferation a substance found in the blood. In other words, while 
the proliferation is a local mechanism, the factor leading to competition is a 
global one . In ecological systems the global character of the limiting factor 
may be related with the fact that the predators move much faster than their 
pray and that their population is proportional to the total prey population. 
Therefore their population constitutes an indirect interaction between distant 
prey locations. We will show that in such cases the effects of inhomogeneity 
are even more important than in the previous sections. As an aside let us 
note that this kind of reactions although formally local (A lion can eat a pray 
only if they meet) are very effective in inducing top down effects; that is the 
regulation of the local pray population according is enforced by their total 
population. 

The overall homogeneity and uniformity assumption used in ODE has 
the benefit of merging into one formalism a large number of biologically 
different homogeneous systems. This advantage turns into a disadvantage 
when multiple spatial scales are involved in the system, as will be described 
in this section. MS on the other hand can describe in a natural way only local 
interactions. It is non-trivial to describe with MS factors acting on widely 
different scales. The solution in this case would be to create a new type of 
HM, which will allow in particular to express the interaction between agents 
and their averages over various ranges. 

A simple example for such system is obtained by modifying the regulation 
mechanism of the system presented in the previous section (Figure 7). There 
the 2 mechanisms proliferation and competition were local interactions. We 
can now replace the local competition by a global competition. Even-though 
from the microscopic point of view the non-locality might look unnatural this 
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is the only way that nonlinearities can arise at the macroscopic level. Indeed 
A term of the form < A > 2 can appear only as the average over terms of the 
type Ai* < A >. Therefore the modification of Eq. |6] would be 

Ai = d 1 S i A i - d 2 Ai - XAi <A t > (9) 

That is we replaced the local competition in the system described in section 
4.1 by a global competition term. Once again we find that passing from local 
to global reactions has crucial implications: By comparing the behavior of 
the system Eq. |9] with the system Eq. [8] one find that they are dramatically 
different in all the parameter regions. Even in the region in which the system 
Eq. |8| is described correctly by an ODE, the ODE are still incapable to 
describe properly the system Eq. ||. Indeed in the | system the entire space 
is empty of As except for a single island centered around the point in which 
the S concentration is maximal (Figure 8). 

We can understand this system by simplifying it to :[] 

A i = A i S i -A l <A I >=A i {S i -<A i >) (10) 

One can see that in every point in which Si is lower than S max the Ai 
density will vanish, while around the point i max in which Si = S max a region 
of high Ai density will appear. Suppose that one starts with all the AiS very 
small; As long as the average concentration of the A's (< Ai >) is lower 
than S max , the local concentration Ai at the point i max will rise (acording to 
Eq. 10). This in turn will raise the value of < Ai >. Eventually < Ai > 



will reach the value of S max . What will be the equilibrium value of the 
AiS in every other point? In every point where Si < S max , one will have 
Ai(S— < Ai >) = Ai(Si — S max ) < 0, which according to [L(] means that 
the concentration at those points will decay to 0. It is clear therefore that 
the entire A population responsible for < A { >= S max is due to a very dense 
island situated around i max - The equality < A { >= S max is not affected 
by the diffusion, which only has the effect of spreading the distribution of 
the A^s in a limited region around i ma x- One sees that the dynamics is 
dominated in this case by the most extreme stochastic local fluctuation in 
the Si distribution, and not by the global or average properties of the S^s. 



1 We rescaled the equation to di = 1, A = 1 and added the death term d-i into the S 
density. This can be done with no loss of generality 
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4.3 Dynamical Ss 

In the previous section we considered a static Si distribution. The intrinsic 
microscopic inhomogeneity of the Ss produced large-scale inhomogeneity in 
the concentration of the As. We will show that these effects survive the in- 
troduction of S diffusion, as long as the maximal proliferation rate of the A's 
is not lower than the diffusion rate of the Ss. If the signaling agent diffusion 
rate of the S is very high, the As will "see" a blurred S concentration and will 
react only to the average S concentration, as in the ODE limit. If however 
the diffusion rate of the S agents is not much higher than the As proliferation 
rate then the As do react to the local concentration of the Ss. The As pro- 
duce peaks in their concentration at the current locations of the S maxima. 
The locations of the S maxima change since the Ss diffuse. Consequently 
the regions of high A concentration will follow them. This will lead to long 
term large (intermittent) fluctuations in the total A concentration (Figure 

9.)- 

In this system the difference between the ODE and MS is not due to 
any imposed inhomogeneity in the starting conditions. It is only due to the 
discreteness of the Ss. 

4.4 Emergence of Complexity 

In the first chapters it has been shown that the discrete character of the 
elementary agents S leads to unexpectedly rich and complex behavior of the 
As population. The complex spatio-temporal emerging structure of the A 
population cannot be predicted by an ODE approach. Moreover in this ap- 
proach one cannot even ask the appropriate questions about the system's 
dynamics. However a system which has only proliferation and death (linear) 
terms leads eventually to a divergent A population. In the present chapter 
we introduced a few types of competition terms. Local competition limits 
the size and maximum concentration of each macroscopic A island. How- 
ever rather than interfering with the emergence of macroscopic features, the 
competition terms turned out to enhance the local and complex character of 
the A islands. Global competition can endow these macroscopic islands with 
an emerging social (or sometime anti-social) behavior. The evolving macro- 
scopic objects may now have apparent "goals": seeking food, multiplying, 
dying or destroying the competing islands. The lifespan of these emerging 
islands is many order of magnitudes longer than the lifespan of their compo- 
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nents (which of course remain non adaptive and totally mechanical during 
the entire history of the system). 

5 Negative Feedback Loops 
5.1 The Survival of the Weakest 

The complex behavior observed in the preceding sections is the result of 
positive feedback loops (autocatalyis). One could think that the difference 
between the ODE and MS descriptions is specific to systems with autocat- 
alytic components. We will show in this section that this is not the case. In 
fact we consider a simple negative feedback situation: The response of the 
immune system to a pathogen invasion. The high pathogen concentration 
leads to an increase in the immune cells concentration which in turn leads to 
a decrease in the pathogen concentration. 

One can express the dynamics of this system using ODEs, in which x is 
the concentration of the pathogen and y is the concentration of the immune 
cells. 

x = —axy (11) 
y = bxy — cy = (bx — c)y (12) 

Eq. [Tl] expresses the fact that upon meeting an immune cell the pathogen is 
destroyed. Eq. [12] expresses the fact that the immune cells proliferate when 
they meet a pathogen with a probability (b), and die with a constant rate 
(c). The solution of the ODE consist in a fast rise of the immune cell 

population, and following it a decrease in the pathogen population, until the 
pathogen population is lower than 0sciU c ationsh - After that point the immune 
cells population starts to decrease. The dynamics of the MS and the ODE 
are the same in this case, with a very notable difference in the final steady 
state of the system: 

• Using ODE one gets 2 types of asymptotic steady states depending on 
the initial condition (x(0),y(0)). More precisely if be* (ln(bc) — 1) — 
ay(0) + bx(0) — c * ln(bx(0)) < both the immune response and the 
pathogen decreases asymptotically to zero. Otherwise the immune re- 
sponse asymptotically decreases to zero, while the pathogen population 
approaches a non-zero asymptotic value. This can be interpreted as a 
constant latent infection. 



5 NEGATIVE FEEDBACK LOOPS 



15 



• In the MS various spatial regions may converge either of the 2 steady 
states described above within a finite time. The regions containing a 
non zero concentration of pathogens (and no immune cells) constitute 
protected reservoirs of pathogens. 

This can explain for example the fact that even people receiving po- 
tent anti HIV drugs for which the average concentration of HIV is 
undetectable for a long period reproduces the disease spontaneously 
when the drugs are stopped. 



5.2 MS Eliminates Spurious Oscillations 

It has been shown in the previous sections that ODE fail to represent emer- 
gent features of dynamical systems, but ODEs can also generate artifacts 
that do not appear in real systems. Such artifacts are exemplified by the fol- 
lowing system consisting 3 cell types: An external antigen, idiotypic (Id) and 
suppressor cells (A-Id). The antigens stimulate the Idiotypic cells. The idio- 
typic cells in turn stimulate the anti idiotypic cells, which destroy idiotypic 
cells. We also introduce a saturation term that regulates the anti idiotypic 
cell population (Figure 12). The ODE describing this system is : 

x = ax — bxy (13) 
y = cxy — dy (14) 

where x is the concentration of the Id cells and y is the concentration of the 
a-Id cells. Such equations have an oscillatory solution around x = -,y = | 
with a period of yda (Fgigure 13). 

The same reactions simulated at the microscopic level do not produce 
oscillations [25]. The MS produces regions containing a high concentration 



of Id cell with very low concentration of A-Id cells around them, or regions 
of A-Id cells that completely destroyed the Id cells around them dominate 
the concentration. The disappearance of the oscillations in the MS is due 
to the de-phasing of the different regions in space. Each regions variate 
independently so that the average yields a constant concentration. When 
one cluster has a maximum concentration of Id cells, other clusters can have 
a low concentration of Id cells. Therefore the oscillatory character of the 
ODE solution is related with the unrealistic assumption that the system is 
precisely homogeneous in space. 
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If the concentration of the Id and Anti Id are low, one can still have in the 
MS system quite significant fluctuations around the average due to the finite 
size of the system relative to the clusters size (Figure 14) (i.e the fact that 
the system contain at every moment a small number of clusters, so that the 
variation of every cluster influences the average behavior). This fluctuations 
involve the total disappearance of agents in certain regions of space. In such 
a situation it might take a macroscopic amount of time until this region is 
"colonized" again by agents immigrating from other clusters. 

This simple system shows that the appearance of oscillations cannot oc- 
cur spontaneously, unless there is an explicit mechanism phasing the whole 
system. This mechanism can be a very high diffusion and/or an external 
input forcing a unified timing over the whole system. Oscillations appear in 
many descriptions of systems using ODE [16|, [17J . This simple system shows 
that these oscillations are many times only an artifact of the model continuity 
assumption, and not a realistic feature of the system. 

6 Discussion 

Dynamic Biological systems were traditionally described using ODE. ODE 
were invented and proved useful in simple physical systems that fullfiled 
certain restrictions: 

• A very large density of interacting particles at every point in space. 

• A high reaction rate. 

• A diffusion rate higher than the local interaction rate. 

• An immediate result of each interaction, i.e. no delay between the 
interaction and its result. 

All of these assumptions are usually false in biological systems. 

• The density of interacting particles is small. For example the number 
of interacting cells at every point is of order of tens. 

• The reaction rate is low - The interaction between cells is a long process 
involving a long chain of molecular sub processes. 
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• The diffusion rate can in many cases be low. The motion of cells cannot 
be described as a diffusive process. The motion is best described as a 
slow directional motion directed by many chemical signals. 

• Most biological processes require the production of new molecules or 
the division of cells. Each of these process can take hours - days. Thus 
there is a delay between an action and its result. 

Moreover the dynamics of physical systems is dominated and in fact de- 
termined by conservation laws. In biology none of the classical quantities 
(mass, number, energymomentum....) is conserved. 

In this work we showed that even in simple cases ODEs fail to describe 
appropriately the dynamics. This is as well: Systems which were too sim- 
ple to present complex emergent features in the framework of ODE turned 
out to present when studied appropriately a wide range of experimentally 
documented features . We used in addition to the straight-forward MS in- 
termediate techniques (HM) that allow to take explicitly into account the 
intermediate relevant dynamical scales between the MS and ODEs. 



6.1 Modeling Autocatalysis 

One of the main sources of complexity is the combination of microscopic in- 
homogeneity and autocatalysis ||27|| . In biological systems there is an inherent 
source of microscopic inhomogeneity- the discreetness of the entities compos- 
ing biological systems. Moreover most biological systems involve prolifera- 
tion which is the basic autocatalitic mechanism. In reality most situations 
are mixed containing both autocatalytic and regulatory mechanisms and can 
be modeled using custom made models describing globally the non catalytic 
parts and explicitly the microscopic details of the autocatlytic parts. 



6.2 Inter-Scale Information Flow 

Autocatalysis results in the flow of information from small scale to large-scale 
feature |28|; from the microscopic details of the system to its macroscopic 



behavior. The rapid evolution of small-scale features to a macroscopic scale 
affects the distribution of the macroscopic variables. In particular the macro- 
scopic variables which are dominated by the sum of a small number of large 
values (instead of on the sum of a large number of small values) have a very 
specific (fractal-power law) space-time behavior. Such effects are enhanced 
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when there is a direct effect of the macroscopic values on the microscopic 
mechanisms. Such a dialog between different scales is the obvious reality in 
biological systems [^| . 

Many biological systems evolve at time scales much longer than their fun- 
damental microscopic time scales. In ODE the emergence of long time scales 
cannot occur naturally, while in MS they can emerge spontaneously. In MS, 
the agents self-organize |30| in collective structures with emergent adaptive 
properties and with dynamical space-time scales of their own. The effective 
dynamics of these emerging objects is much slower than the one of their com- 
ponents (and this fact constitutes their very definition Analyzing the 
system at the level of the emerging objects require automatically the tran- 
sition to longer time scales and larger spatial scale. Each level corresponds 
to a new space-time scale in a hierarchy of complementary representations 
of the same complex macroscopic system. 

We found in this paper that the core mechanism for the emergence of 
collective macroscopic objects is the interplay between autocatalysis and the 
discrete microscopic structure of the components of the system. This mecha- 
nism is totally missed by the mean field description, which predicts uniform 
extinction in conditions in which the real system emerges actually very lively 
collective adaptive objects. 
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Figure 1: The reactions taking place in a system of proliferating cells. 1) 
Whenever an A agent meets an S agent, a new A agent is created. A + S — > 
S + A + A, with a rate of a?i 2) Each A agent has a constant probability <i 2 of 
dying. A — > $. The Ordinary Differential Equation (ODE) approximation 
in the continuum limit for these 2 reactions is: A = (< S > *d\ — d-z) * A. 
In the left figure full arrow represent activation, and dotted arrows represent 
destruction. 
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Figure 2: The solutions of equation 3. An exponentially decaying solution 
for (< S > *di — ri 2 ) < , an exponentially growing solution for (< S > 
*di — d 2 ) > and a constant solution for (< S > *di — d 2 ) = 0. 
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concentration of S agent. 



Log concentration of A agent. 
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Figure 3: Comparison between the MS and ODE description of the system 
described in figure 1. The A drawing is the local concentration of S agents 
at every latice site. The B drawing represents the log of the A agent concen- 
tration at every site. We use a log representation in order to show the large 
variance in the local concentration of the A agents. The C drawing shows 
the average A agent concentration in the MS (interrupted lines) and in the 
differential equation model (circles). 
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Log A agent, concentration in MS 



Log A agent, concentration in HM 
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Figure 4: A comparison between the concentration of A agent in a MS and 
in HM model. One can see that the concentration is similar in the 2 models, 
if one uses appropriate parameters. Comparison between a HM, an MS and 
an ODE. In the lower drawing we show that hybrid models can produce 
both the results of the microscopic representation, and the results of the 
differential equations. Figure A is the comparison between ODE and a HM 
for an homogeneous S concentration. Figure 4B is the comparison between 
an MS and a HM for a non- homogeneous ^concentration. 
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Figure 5: Reaction scheme for a system containing local competition. This 
system has constant proliferation and death rates as in Figure 1. However it 
has an extra mechanism of local competition: Whenever two A agents meet 
there is a probability that one of them will be destroyed A + A— > A 
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concentration of S agent. Log concentration of A agent. 
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Figure 6: Comparison between ODE and MS for the system defined by figure 
5. This system has an average negative growth rate yielding an exponential 
decay in the ODE description as in the previous section. In the MS we 
observe a limited growth of the A islands, which is due to the competition 
between As. Thus even regions in which the local A concentration would 
grow are limited to a finite concentration. 
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Figure 7: Global competition reaction scheme. The mechanisms of prolif- 
eration and death are similar in this system to the one described in figure 1. 
The extra mechanism in this system is the global competition: Each A agent 
has an extra death rate proportional to the average concentration of the As. 
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Figure 8: Comparison between ODE and HM for a system with an average 
negative growth rate and a global competition term (as defined in figure 7). 
In this model a single large island with a very high A concentration appears. 
The competition with this island destroys all the other small A islands. The 
A and B drawings show the log of the A agents concentration at different 
times. The C drawing shows the comparison between the ODE and the HM. 
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Log concentration of A agent at time 20 
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Figure 9: The model defined in figure 7 in a regime with high S diffusion 
rate. When the S diffusion rate is very high, a few large A islands are 
created. As opossed to figure 8, the largest A island does not dominate. The 
various islands emerge an adaptive behavior: They look for the maxima of 
the S concentration, split, merge and die. The top windows shows different 
snapshots of the A distribution. The small islands are destroyed while the 
large islands grow. The bottom window is the a long term evolution of this 
system showing intermitent fluctuations. 
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Figure 10: The reactions scheme of a system describing the destruction of a 
pathogen by the immune system. 
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Figure 11: The simulation by a MS of a system describing the down regula- 
tion of pathogen by immune cells. In contradiction with the ODE description 
that predicts an uniform steady state, the MS predicts the formation of pro- 
tected zones in which the pathogen can survive. The C and D drawings 
represents 2 cases predicted by the ODE description: In the C drawing the 
pathogen is not destroyed, and in the D drawing the immune system manages 
to destroy the pathogen. The A drawing describes the pathogen distribution 
simluated by the MS. One can see it contains protected regions containing 
pathogens surrounded by empty regions. 
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Figure 12: Reaction scheme for a system describing the regulation of idio- 
typic cells by anti-idiotypic cells. 
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Figure 13: A comparison between the MS and the ODE for a system contain- 
ing idiotypic and anti idiotypic cells. The ODE leads to oscilations. These 
oscialtions disapear in the MS when the random reactions are taken into 
acount. The MS predicts that the Id and A-Id cells concentration decays to 
a constant. 
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Figure 14: Large fluctuations can appear in the MS description of a system 
with Id and anti Id cells. Most of the contribution for the concentration is 
in very localized islands of very high concentration. These leads to the very 
sharp fluctuations. This system reaches a final steady state in which space 
is occupied by large disjoint regions of Id and respectively A-Id cells. Each 
such region has non stationary dynamics. 



